HIGH-ORDER FINITE VOLUME SCHEMES FOR LAYERED 
ATMOSPHERIC MODELS 
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Abstract. We present a numerical scheme for the solution of a class of atmospheric mod- 
els where high horizontal resolution is required while a coarser vertical structure is allowed. The 
proposed scheme considers a layering procedure for the original set of equations, and the use of 
high-order ADER finite volume schemes for the solution of the system of balance laws arising from 
the dimensional reduction procedure. We present several types of layering based upon Galerkin 
discretizations of the vertical structure, and we study the effect of incrementing the order of hori- 
zontal approximation. Numerical experiments for the computational validation of the convergence 
of the scheme together with the study of physical phenomena are performed over 2D linear advective 
models, including a set of equations for an isothermal atmosphere. 
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1. Introduction. When we study phenomena in the atmosphere, they are often 
confined to certain layers, for example the atmospheric boundary layer (ABL) and the 
stratosphere. Models of atmospheric phenomena often reflects this layering. Examples 
of use of layered models are: air pollution models (see for example the early study 
in [T]), moisture and precipitation studies [3], hurricane modelling [TJ, large scale 
atmospheric dynamics and climate models. 

There are differences between a conventional vertical discretization of a 3D atmo- 
spheric model, and a layered model. The most important difference is that layered 
models are focused on specific properties of certain layers, whereas the purpose of a 
conventional vertical discretization is to express numerically the vertical variation of 
the complete atmosphere. Vertical discretization do treat say the ABL specifically, 
but only by having sufficiently man vertical layers so as to resolve the phenomena of 
interest. 

The numerical schemes often used in atmospheric models are some form of conser- 
vative difference schemes in the vertical and second order centered difference schemes 
on a staggered grid in the horizontal. The global weather model at the European 
Centre for Medium Range Forecasts (ECMWF) introduced a finite element scheme in 
the vertical relatively recently, and finite volume schemes has also been constructed 

m- 

What we want to communicate here is an investigation of using of a high- order 
(> 2) finite volume discretization in the horizontal combined with several different 
schemes in the vertical. The purpose is to find out how efficient and accurate such 
discretizations are. The horizontal discretization that we apply is a finite- volume high- 
order Godunov scheme called ADER (Arbitrary high-order DERivatives) , proposed 
by Toro and Titarev (see [2^ and [24], among others). The ADER scheme is actually 
a family of schemes with increasing order. Hence there is a trade-off between order 
and computational complexity. The schemes have been shown to work well in a series 
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of applications, see for example [5] and [37]; related schemes, applied to shallow water 
models, are presented in [i], [B] and [S]. 

A layered model is obtained by some form of averaging of the governing equations 
in the specified layer. This procedure then leads to a set of equations in the horizontal, 
for each layer; the equations in the layers are then coupled. An example of such 
procedure for the linearized 3D Euler equations with an orthogonal set of finite element 
hat functions can be seen in obtaining what is called a "2.5D" or "1.5D" model. 
We extend this approach by not only considering piecewise linear finite elements, but 
also a pseudosprectal method, a layering based on discontinuous Galerkin methods 
and a method based on vertical averaging imposing conservation properties. Such 
investigation in the layering procedure arises from the fact that apart from the physical 
basis of a layering, there is a potential saving in computational time, since only a few 
layers are required and application of numerical schemes specifically aimed at 2D or 
ID problems can be used. 

One of the motivations for our work is its application to data assimilation, more 
specifically ensemble- type of methods as the Ensemble Kalman filter (EnKF) . For an 
overview, see e.g. [I5]- The EnKF procedure involves an ensemble of model runs 
between each observation time point, and preferable the size of this ensemble should 
be as large as practically possible. This implies that model runs should be as efficient 
as possible, and thus, the use of dimensionally reduced models is highly desirable. It 
is worth to emphasize that our focus on data assimilation is somewhat different from 
conventional data assimilation in weather models ( treated in e.g. [TH] and [5]). We 
are concentrating on high resolution models where there are a few observations in the 
ABL, and we want to study the effect of these observations in the other layers of the 
atmosphere. A typical application where such a situation is relevant, is windpower 
modelling. See [21] for a detailed description. 

To computationally study convergence, and to assess the performance of the 
method, we first use a linear 2D advection where we investigate the properties of 
the numerical schemes. Then we proceed to approximate a linear, non-hydrostatic 
atmospheric model with the same family of schemes; apart from convergence issues, 
we present a physical test related to the propagation of inertia-gravity waves of non- 
hydrostatic scale in a periodic channel. 

This paper is organized as follows: In section 2 we describe in detail the verical 
discretization (layering) methods for a 2D linear advection-reaction problem. Section 
3 contains a brief description of the ADER scheme for horizontal discretization. Nu- 
merical tests validating convergence rates and other relevant physical phenomena for 
2D linear advection and for a 2D linear atmospheric model are shown. Outlook and 
summary is included in section 5. 

2. Vertical discretization. We consider a linear, constant-coefficient, two- 
dimensional advection-reaction model given by 

dtq + adxq + bdzq + sq — 0. (2.1) 

The problem is spatially defined over a rectangular domain fl — {{x,z) G [0,ia;] x 
[0,i?2]} together with compatible boundary and initial conditions. We study four dif- 
ferent methods to reduce the dimension of the model by discretizing the vertical 
component. It is relevant to note at this point that the inclusion of an additional hor- 
izontal dimension poses no additional difficulties, while the procedure strongly relies 
on the linearity of the vertical structure. The first method that we discuss is a vertical 
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average approach based on vertically constant levels, defining the interface interac- 
tion in a way that allows us to enforce some conservation-like properties. Secondly, 
we approximate using a continuous finite element expansion, and in a familiar 
formulation we include layering via a pseudospectral collocation method. We finally 
remove the continuity assumption in the finite element setting and we perform lay- 
ering with piecewise P*^ and and P"'^ space-discontinuous Galerkin methods; note that 
in this latter case, the technique is applicable to nonlinear problems in a tractable 
way. For more details on the formulation and application of such techniques for the 
numerical approximation of time-dependent problems please refer to |25l and [3] (and 
references therein). 

2.1. A "conservative" vertical averaging (CVA) approach . We start by 
defining a set of uniformly distributed vertical levels Zi^ i — . . . N , with zq — and 
zn = Hz ■ For every interval h = [zi-i, Zi] we define the vertically averaged quantity 

9« = r / q{x,z,t)dz, (2.2) 



where h = H^/N . We consider a local equation for every and we proceed by 
vertically integrating it, obtaining an exact formula for q^. 

dtQi + adxQi + ^b{ q{zi) - q{zi^i) } + s = 0, (2.3) 

We focus our attention on the approximation of the interface values q{zi) and q{zi^i). 
We proceed by assuming a linear vertical variation between these two quantities within 
li and imposing a consistency condition with the vertical average 

-J (z - z.,) + q{zi^i)dz = q^, (2.4) 

which after integration is reduced to 

q{zi) + q{z,-i) K 2q,. (2.5) 

This way to approximate interface values gives a coupling between the N levels, and 
therefore the original 2D system is reduced to a one-dimensional system with a source 
term that mimics the vertical interaction; for instance, if we impose q{zo) = in order 
to close the set of equations, this approach leads to 

dtqk + ad^qk+bCqk + sqk =0, Qk = [qo ■ ■ ■ QnY, (2.6) 

with 

2 if I = J, 

C^J = f^ (-1)^+^4 if*>j, (2.7) 

if i < j. 



Note that such formulation generates a source term C which, for a fixed mesh pa- 
rameter h, has only positive eigenvalues A = ^. 
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2.2. Piecewise linear finite element (FEM) expansion. A different layered 
formulation can be obtained via a Galerkin ansatz in the vertical direction. Proceed- 
ing in a similar manner as in the method of lines, we assume that q{x,z,t) can be 
approximated by an expansion of basis functions in the vertical direction with space- 
time dependent coefficients, 

g(a;,z,t) w g,, = ^gfe(x,0*fcW- (2.8) 

k 

For a fixed pair (x, t), we denote by i^(0; Hz) the space of square integrable functions 
in z with its associated inner product 

(u,«)l2(o;//,) = / u{x,z,t)v{x,z,t)dz, (2.9) 
Jo 

and standard space H^{0;Hz) is defined upon the inner product 

{u,v)m{0;H,) = {u,v)l^O;H,) + {d^U, 5^ W ) £2 (0;//, ) • (2-10) 

We uniformly partition the interval [0, H^] into N elements li ~ [zi-i, Zi] i ~ 1 . . . N 
and mesh parameter h = Hz/N. The approximation space Vh C H^(0; Hz) is defined 
as 

V,, = {ve H\0; Hz) : v\u e F\ i = I, . . . , N}. (2.11) 

The set { $^(0) '^^ ^at functions constitutes a basis of the finite dimensional space 
Vh ■ We write a weak formulation for (|2.ip , 
Find qh & Vh such that 

{{dt + adx + s)qh,v)L2(^Q.H,) + {bdzqh,v)L2(o;H,) = , Vw G V/^. (2.12) 

In the same manner as in our first approach, this formulation is formally reduced to 
a ID linear hyperbolic system for the expansion coefficients {qk}'- 

M{dt qk+adxqk + sqk) + bKqk^O, qk ^ [qo ■ ■ ■ QnY , (2.13) 

where 

M,,j- = ($,,$j )i2(o.^^), K,^j ^ {^z<^>^,<^>J)L^o■M,)■ (2-14) 

Since the inverse of M has to be computed only once, the simplified version reads 

dtqk+ad^qk + bCqk+ sqk = 0, C = M'^K. (2.15) 

Although the two layering approaches so far presented lead to one-dimensional 
hyperbolic systems of balance laws with a similar structure, the representation of an 
approximate two-dimensional solution is radically different. The first method will 
allow to recover piecewise constant average values for each vertical level, while the 
finite element method, with this particular choice of basis functions will be able to 
recover a piecewise linear and continuous approximate solution of the two-dimensional 
problem. In the same spirit, the forthcoming section develops an expansion in the z 
direction such that a polynomial of arbitrary degree can be recovered. 
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2.3. Pseudospectral collocation (PS) method. The formulation of a pseu- 
dospectral method can be interpreted as a Galerkin approximation when a weight in 
the inner product is introduced and specific quadrature points are chosen. We be- 
gin by mapping the domain il into ftp = [0, Lx\ x [—1,1], where we consider again an 
expansion of the form Qp = J2k Ik i^^ {z) , with a basis given by the set of Tscheby- 
chev polynomials {Ti~{z)} up to degree p; it is well-known that such polynomials are 
an orthogonal set in under the weighted inner product given by 



(w, u) = / ■ 7((t- 7- 7)(7- z dy. (2.16) 

J-1 Vl — z'^ 

Given p, we define p + 1 collocation points Zi which in our case are optimally chosen 
as the zeroes of the Tschebychev polynomial Tp+i(z), namely 



- cos 



{2i + l)7r 
2(P+1) 



z = 0,...,p. (2.17) 



We obtain a set of p + 1 equations by imposing an exact solution at every collocation 
point, i.e., 

p 

{dt + adx + bd^ + s)"^ qk{x,t)^{z) ^ , at every z = z^. (2-18) 

The reduced system now is given by 

dtqk + adxqk + bCqk + SQk ^0, qk = [qo ■ ■ ■ qpf , (2-19) 

with 

M,j =T,(z,), = 9,T, (z,), and C = M^^K. (2.20) 

2.4. Discontinuous Galerkin (DGPO - DGPl) approximation. Proceed- 
ing in a similar way as in the piecewise continuous and linear finite element approach, 
after partitioning the vertical domain we consider a local problem in the interval 
li ~ [zi-i, Zi], which is mapped into the reference element Iref = [—1,1] via 

z='-^+4 (2.21) 

In such domain we consider a local expansion 

g,(a;,z,t) «^c,,,;(x,t)$,(z), (2.22) 
i 

where the basis { $j(z) }^^q is chosen as the set of Legendre polynomials up to degree 
p which is orthogonal over [—1,1]. If p = 1, then qi{x,z,t) — co,i +^ci^i and after 
multiplication by test functions in the space of polynomials of degree 1 over [—1,1], 
P} ^ , and integration by parts, we write the following weak formulation (with the 
same product as in (|2.9I) but in the interval [— 1; 1] ) for the local problem : Find 
qi e P}^_^^ sucii that 

{{dt +adx + s) gi,w)L2(_i.i) +vhq\\^- {bqi, 9^ =0 Vw G 

(2.23) 
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Note that the basis in the expansion are no longer of compact support as in the 
piecewise hnear continuous finite element approach, and therefore there is no conti- 
nuity enforcement; every local problem is treated separately. However, we do enforce 
coupling between the different layers by a suitable approximation of the interface 
values arising in the term vbqi\l^j^; in this article we approximate this interface flux 
via upwind considerations. Before focusing our attention on this issue, we recall that 
testing is done with v = 1 and v = leading to the following set of equations for the 
coefficients co,i and ci^i at every element: 

dt co,i + adx co,i + sco,i + ■^{bq{z ^ Zi) - bq{z = Zi-i) } = 0, (2.24) 
3 

dt ci,i + adx ci,i + sci^i + - {bq{z = z,;) + bq{z = Zi-i) - 2bcoA } = 0. (2.25) 

This system of equations can be closed by suitable boundary conditions on vertical 
extremes of the domain; the remaining issue is henceforth the approximation of the 
interfaces values. We use an upwind flux, that for 6 > leads to: 

bq{z = z,) = bq,{i = l) = fe(co,. + ci^,) (2.26) 



bq{z = Zi_i) = bqi^i{^ = 1) = &(co,j-i +ci,i-i). (2.27) 

As in all our previous approaches, the dimensionally reduced set of equations is a 
one-dimensional system of hyperbolic balance laws, consisting in this case in two 
equations determining the local coefficients cp^i and ci^i for every element; there is an 
increasing computational cost depending on the local approximation degree compared 
with continuous formulations. 

The next section takes as a starting point the resolution of a hyperbolic system 
of balance laws in the form of those obtained by the layering procedures so far in- 
troduced; they all share a single horizontal advection speed while a stable reactive 
source establishes coupling and transmission between the different vertical levels. The 
same approach can be extended to systems where different horizontal advection speeds 
are present, as we will show in the examples, but since it is possible to work in the 
characteristic variable field they key points remain in the resolution of this simpler 
problem. 

3. The ADER approach for high horizontal resolution. This section con- 
sists of a concise version of the original scheme formulation; it makes extensive use of 
the linearity of the system while the general formulation has been developed by con- 
servation laws in a wide framework. For more details about the general formulation 
please refer to [26l EH EH] and references therein. We consider a system of balance 
laws in one spatial dimension given by 

dtq + Adxq + Sq^O, (3.1) 

along with initial and boundary conditions. We will assume A diagonal, and S 
a definite semipositive matrix. We partition the horizontal domain into intervals 
Ej = [xj_if2i 2;j+i/2], and we consider a control volume in the x—t space Ej x [t", 
with dimensions Ax — Xjj^xj2 ~ Xj^xjii At = — t". Intcgrating p.ip over the 
control volume leads to 

-1?-^ - - (3.2) 
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where 



1 , . . 11 



/+1/2 = ^ y ^ g(a;j+i/2,T")dT, ^1 ^ '^'^ J J Sq{x,T)dxdT, (3.3) 
which provides a marching exact formula for the spatially averaged quantity 



1 r^]+i/2 

q] = -^ q{x,e)dx. (3.4) 



The ADER scheme is a numerical method that generates a high-order approximation 
of the two integrals appearing in (I3.3P and consequently, a high-order approximation 
of the averaged state value. The key ingredients of this scheme applied to this system 
are: 1) a reconstruction procedure that allows us to recover high-order degree poly- 
nomials at every element given the set of averaged values, for 2) the formulation of 
a generalized Riemann problem for, 3) the generation of Taylor series expansions for 
the approximation of the required states. After integration is performed, the resulting 
scheme is of order m, both in space and time, where m — 1 corresponds to the polyno- 
mial degree used in the reconstruction step, provided that the integration procedures 
are of at least the same order of accuracy; in our case, time integrals will be treated 
exactly, while the space integral remaining in the source term will be approximated by 
a five point Gaussian quadrature (in this article we study schemes up to order 4) . For 
the reconstruction step, we the use WENO procedure proposed in TU], which slightly 
differs from the original version |221 117) , as it recovers a single polynomial for every 
cell, which can be used either for extrapolation at the boundary or evaluation inside 
the cell. As we will be dealing with systems of balance laws, this latter characteristic 
is useful for the approximation of the source term. Henceforth we denote by ADERm 
to the scheme that uses a to — 1 polynomial degree reconstruction. Wc now focus 
on the task of approximation of the fluxes in (j3.2p , which requires the approximation 
of the state q at the interface points in every horizontal element, i.e, q{xj^i^2T''') 
and q{xj-i/2,T). After every time step inicialization at i", being performed the re- 
construction step, for every element we dispose of a reconstruction polynomial Pj{x) 
of degree to — 1. At every cell interface we pose the following generalized Riemann 
problem: 

PDE: dtq+ Ad^q + Sq^Q, (3.5) 

ic:qix.,)^{ If^r;^^:^^^. IIT'^ (3.6) 

The ADER approach links the resolution of this generalized Riemann problem with 
the approximation of the interface state qix^j^xji^T). Assuming t sufficiently small, 
we write a Taylor expansion in time for the interface state 

m— 1 ^ 

'7(2^,+i/2,t) «g(a.j+i/2,0+)+ ^[atWg(a;,+i/2,0+)]^, (3.7) 

k=\ 

where the time derivatives are replaced with space derivatives via the Cauchy-Kowalevski 
procedure ( also know as Lax-Wendroff procedure [20) ) procedure. In the simplest 
case, if ^ is a diagonal matrix with a single element, the k — th time derivative holds 
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the binomial formula 

dt^''\{x,t)^{-l)\Ad., + Sfq{x,t), (3.8) 



i=0 



{-lfJ2(jA'S''-'dJ''^q, fc = 0,...,m-l. (3.9) 



The interface state and its spatial derivatives are approximated by the solution at 
(x — Xj_^_i/2)/T—0 of the following set of classical Riemann problems derived from 
(|531) : 

PDF : dt 8:,^''^+ Ad,,^''\ = 0, fc = 0,...,m-l (3.10) 

IC ■ d ^^^q(x 0) = I '•'"^ j+1/2) = ^^h]{xj+i,2): X < .Tj+i/a 

I d^'^''''qRiXj + i/2) ^ d^^^'>Pj + i{Xj + i/2), X>Xj + i/2 

(3.11) 

There are two things worth to point out at this stage: First, that given A a diagonal 
matrix with a single entry, the formula p. 81) uses the fact that A and S commute; 
the case with multiple entries in the diagonal is still tractable even though it is not 
possible to make use of the binomial formula; this same diagonal structure allows 
an exact solution of p.lOp component-wise where it is reduced to a linear advection 
problem, while the more general formulation of the ADER scheme include the use 
of approximate Riemann solvers at this point. Secondly, we note that although the 
set of classical Riemann problems used to approximate q{xj^i/2,0~^) and its space 
derivatives neglect the presence of the source term, it is still captured in the final 
Taylor expansion 



q{x,+,/2,T) « g(a;,+i/2,0+) + E 1 E J -^''5'"'^- ('='g(a;,+i/2, 0+) 

k=l U=0 ^ ^ 

We insert this formula into the flux p.3p and integrate, yielding 

-4 Lx^.^,/,,0+)+E7E(J)-4'5'=-9.Wg(x,+V2,0+)| 



fc! ■ 
(3.12) 



k=l k i=0 



(3.13) 

The remaining issue on the scheme is the approximation of the source term in- 
tegral Sj] as we have previously seen, since we write Taylor expansions in time and 
the source is linear, we don't have to worry about time integration as it can be done 
exactly; the space integral is approximated by a Gaussian quadrature rule, therefore 



1 1 



r^tr^, +1/2 II ^ 

/ / Sq{x,T)dxdT -— / ^WaS q{Xa,T)dT (3.14) 

Jo ''Xj_i/2 •'O a=l 



where Ng stands for the number of Gauss points {2^0}^=! used in the integration rule; 
its associated weights are denoted by Wa- Within the j — th element, we write a Taylor 
series in time and replace the time derivatives by space derivatives in a similar way as 
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it has been done for the calculation of the fluxes; therefore, inside the element Ej, in 
a given Gauss point Xq. the Taylor formula (I3.12p holds replacing a;j+i/2 by Xa. What 
is different however, is that the state values dx ''^'-'^(xq,, 0+) are not obtained from the 
solution of a Riemann problem, but from the solution of the Cauchy problem at the 
j-th cell 

PDE : dt dx'-''U+ Adx'^''U = 0, A; = 0,...,m-1 (3.15) 

IC : a,«g(a;,0)= dj'''>p,{x). (3.16) 

Performing the whole procedure (reconstruction, flux and source term calculation, 
and update) leads to a scheme of ni—th order of accuracy both in space and time, with 
an optimal CFL condition of CFL < 1 without source term; in certain cases, as we will 
see, the source term does not affect the stability condition and can consequently be 
omitted when determining the time step, while in other situations its incorporation 
restricts the time stepping in a way that other solutions, as splitting schemes are 
required; in any case, the ADER method provides an accurate and efficient way to 
treat the advective character of the studied equations. 

We now perform numerical tests for the global scheme, i.e., the numerical proce- 
dure that includes a low-order layering in the vertical direction generating a system of 
balance laws that is approximated by an ADERm scheme. The purpose of the tests 
are twofold: one on hand they validate the accuracy of the proposed ADERm scheme 
in the presence of different source terms, while they also study the performance of the 
global scheme and its properties. 

4. Numerical experiments. In this section we consider two models in order to 
test the scheme presented in this article. For both cases we consider a two-dimensional 
X — z domain n = {(x, z) S [0, L] x [0, H]}. In every case numerical convergence rates 
are shown in order to assess the accuracy of the scheme, while other relevant tests 
dealing with stability, conservation and physical validation are performed. Most of 
the code concerning the reconstruction procedure, as well as flux and source term 
calculations were coded and vectorized in FDRTRAN95; simulations were performed 
using an Intel Core 2 Duo 2.50 Ghz processor in a computer with 4 Gb RAM. 

4.1. 2D linear advection. The simplest case in which our approach is meaning- 
ful consists on the advection of a quantity q = q{x, z, t) through 51. This phenomenon 
is represented by the linear, two dimensional advection equation 

dtq + Udxq + Wd,q^Q, (4.1) 

where the pair (t/, W) is a constant two-dimensional velocity field. The initial con- 
dition q(x, z,0) is denoted by qo(a;, z). Well-posedness of this problem requires that 
for U and W positive, q is specified at x = and z = 0; we set q{Q, z, t) = q{L, z, t) 
and q{x, 0, t) = qo{x, 0). Given such model, we first perform vertical layering via the 
methods stated in section 2; as it was previously shown, every case leads to a system 

dtqk + Udxqk + Cqu^Q. (4.2) 

where the resulting source term is obtained by replacing h = W . This hyperbolic 
system is solved with the ADER approach described in section 3. 
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4.1.1. Convergence studies. We first perfom numerical tests in order to val- 
idate the convergence rates of the ADER method for a system in the presence of a 
source term. Every source term in this case is stable in the sense that its eigenvalues 
are contained in the left half-plane. We set = Hz = 1, U = W = 0.01 and an 
initial condition 



Time stepping is established as At ~ CFL (Ax/U), with a CFL number chosen to 
be 0.8, being the simulations stable with CFL up to 1; note that the time stepping is 
computed as for an homogeneous linear advection problem in one dimension. Final 
time for simulation is t = 50[s]. For 10 vertical levels (or collocation points in the 
case of the PS source term), and with ADER order varying from 2 to 4 we obtain the 
results shown in Tables 14.1114.21 where relative errors are computed in space at final 
time of simulation with respect to a reference solution obtained with 1024 horizontal 
elements. Predicted convergence rates for smooth solutions, in the presence of every 
source term proposed in section 2, are achieved in both Li and norms. 

4.1.2. Advection of a gaussian initial condition and averaging. In a sec- 
ond numerical experiment we consider the same spatial domain and velocity field, but 
an initial condition 



and a simulation time of 70 [s] (with 100 x 100 elements). We measure the Loo norm 
of the 2D solution in order to discuss the diffusive properties of the global scheme. 
We observe that the application of the ADER method in the horizontal is not causing 
major difussion when orders 3 or 4 are used, while the choice of the source term is 
relevant. Figure HTTl shows that both FEM and PS source terms lead to schemes that 
are not diffusive at all while exhibiting an oscillation due to the transmission of the 
maximum value from one layer to another; the CVA approach shows considerable 
diffusion while the DGPO scheme with upwind flux is the most diffusive one, wich 
was expected since it is the only first order method. Another aspect to be observed 
from these graphs is the ability to capture the correct advection speed; it can be 
noted that the CVA approach overestimates the physical speed, while the remaining 
methods clearly show a decay at t = 50 [s] that corresponds to the expected time for 
the top of the advected initial condition to leave the physical domain. The oscillations 
observed afterwards in the FEM source term case are due to a spurious reflection in 
the upper boundary. 

We conclude this example by analyzing a setting that is relevant in a data as- 
similation context. We consider the same Gaussian initial condition with U = 0.01, 
W — 0.001, but we also consider a subdomain ilg = [0, 1] x [0.5, 0.6] where we vertically 
average the solution. This leads to a quantity 



which is computed upon the approximated solution. Simulation time is now 200 [s] 
with 100 horizontal elements and 20 vertical elements. We show the evolution of Qs, 
which is advected at the same speed, and we show absolute residuals of the averaged 
values obtained from approximated solutions and the exact average value. Figure 1?^ 



q(){x,z) 



sm(107ra;) sin{2nz). 



(4.3) 




(4.4) 




(4.5) 
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Table 4.1 

Linear advection: convergence rates and relative error (% ) for the ADERm—th horizontal order 
scheme with a conservative vertically averaged (ADERm-CVA) and piecewise constant discontinuous 
Garlerkin source terms (ADERm-DGPO) . 



Method N L^o error Loo order Li error Li order 



ADER2-CVA 



ADER3-CVA 



ADER4-CVA 



ADER2-DGP0 



ADER3-DGP0 



ADER4-DGP0 



32 


3 


.3429e- 


•002 






7, 


.5757e- 


•003 






64 


1 


,2924e- 


-002 


1 


.3711 


2, 


.9740c- 


•003 


1, 


.9893 


128 


5, 


.0697e- 


-003 


1 


.3500 


7, 


.59536- 


•004 


1, 


.9692 


256 


2 


.OlOle- 


■003 


1 


.3347 


1, 


.9504e- 


•004 


1, 


.9614 




Q 

o. 






1 

± . 








■uuo 


1 

L . 


Q4SQ 


32 


8, 


.2992e- 


-004 






1, 


.1808c- 


•002 






64 


9 


.7975c- 


-005 


3 


.0825 


6, 


.2395c- 


•005 


3, 


.0906 


128 


1, 


.1854c- 


-005 


3 


.0471 


7, 


.5473c- 


•006 


3, 


.0474 


256 


1 


.4568c- 


-006 


3 


.0244 


9, 


.2749c- 


•007 


3, 


.0246 




1 

L , 






Q 

o. 


ni 9*^ 


1 

J. , 




.nn? 


Oi 


ni 9"^ 


32 


4 


.6056c- 


-005 






3, 


.0473c- 


•005 






64 


2 


.9604c- 


-006 


3 


.9595 


1, 


.8899c- 


•006 


4, 


.0112 


128 


1, 


.8480c- 


-007 


4, 


.0018 


1, 


.1738c- 


•007 


4, 


.0091 


256 


1 


.1457c- 


-008 


4, 


.0117 


7, 


.2826c- 


•009 


4, 


.0106 


512 


6 


.7165e- 


-010 


4, 


.0923 


4, 


.2755e- 


•010 


4, 


.0903 


32 


2, 


.3473e- 


-002 






5, 


.9531e- 


003 






64 


9 


.7303e- 


-003 


1 


.9893 


1, 


.5895e- 


003 


1, 


.9051 


128 


3 


.9960c- 


-003 


1, 


.9692 


4, 


.2553c- 


•004 


1, 


.9012 


256 


1 


.6320e- 


■003 


1 


.9614 


1, 


.1495e- 


■004 


1, 


.8883 


512 


6 


.6328e- 


■004 


1 


.9489 


3 


.0676e- 


■005 


1, 


.9058 


32 


4, 


.2071e- 


■004 






2, 


.6867e- 


■004 






64 


5 


.2219e- 


■005 


3 


.0906 


3 


.3291e- 


■005 


3 


.0126 


128 


6, 


.4997c- 


-006 


3 


.0474 


4, 


.1373c- 


•006 


3, 


.0084 


256 


8, 


.1006e- 


■007 


3 


.0246 


5, 


.1572e- 


■007 


3, 


.0040 


512 


1 


.Ollle- 


■007 


3 


.0123 


6, 


.4369e- 


■008 


3 


.0022 


32 


8, 


.2443e- 


■006 






4, 


.9867e- 


■006 






64 


5 


.7959e- 


■007 


4, 


.0112 


3, 


.6692e- 


■007 


3, 


.7646 


128 


3 


.8606e- 


-008 


4, 


.0091 


2, 


.4559c- 


•008 


3, 


.9011 


256 


2, 


.4817e- 


■009 


4, 


.0106 


1, 


.5797e- 


■009 


3, 


.9586 


512 


1 


.4823e- 


■010 


4, 


.0903 


9, 


.4370e- 


■Oil 


4, 


.0652 
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Table 4.2 

Linear advection: convergence rates and relative error (Yc) for the ADERm — th horizontal 
order scheme with a piecewise linear finite element (ADERm-FEM) and pseudospectral source terms 
(ADERm-PS). 



Method N L^o error L^o order Li error Li order 



ADER2-FEM 



ADER3-FEM 



ADER4-FEM 



ADER2-PS 



ADER3-PS 



ADER4-PS 



32 


2, 


.6751e- 


002 






7 


.5757e- 


-003 






64 


1, 


.1372e- 


■002 


1, 


.2341 


1 


.9927e- 


-003 


1 


.9267 


128 


4, 


.7237e- 


■003 


1, 


.2675 


5 


.39436- 


-004 


1 


.8852 


256 


1, 


.93986- 


■003 


1, 


.2840 


1 


.45906- 


•004 


1 


.8865 




7 






1 

L . 


9047 


■3 

O. 




•uuo 


1 

± . 


QnQ7 
.yuy 1 


32 


7, 


.2227c- 


■004 






4, 


.6239c- 


-004 






64 


9, 


.07430- 


■005 


2, 


.9927 


5, 


.7778c- 


-005 


3 


.0005 


128 


1, 


.13300- 


■005 


3, 


.0017 


7 


.2156c- 


-006 


3 


.0013 


256 


1, 


.41530- 


■006 


3, 


.0009 


9 


.0102c- 


-007 


3, 


.0015 


111 9 


1 

L . 




.nn? 


o. 


nnns 


1 

J. . 




nn7 


Q 

o. 




32 


2, 


.25790- 


■005 






1, 


.4252c- 


-005 






64 


1, 


.45530- 


■006 


3, 


.9556 


9 


.2260c- 


-007 


3 


.9493 


128 


9, 


.27830- 


■008 


3, 


.9713 


5, 


.9028c- 


-008 


3 


.9662 


256 


5, 


.83910- 


■009 


3, 


.9900 


3 


.7170c- 


-009 


3 


.9892 


512 


3, 


.4532e- 


■010 


4, 


.0798 


2, 


.19846- 


-010 


4, 


.0796 


32 


9, 


.3519e- 


003 






4, 


.69196- 


-003 






64 


2, 


.6984e- 


003 


1, 


.7931 


1 


.17436- 


-003 


1 


.9983 


128 


9, 


.64310- 


■004 


1, 


.4846 


2, 


.9331c- 


-004 


2, 


.0014 


256 


3, 


.6181e- 


•004 


1, 


.4143 


7, 


.37976- 


•005 


1 


.9908 


512 


1 


.3932e- 


•004 


1, 


.3768 


1 


.87116- 


•005 


1 


.9796 


32 


1, 


.9786e- 


•004 






1 


.26546- 


•004 






64 


2, 


.3921e- 


•005 


3, 


.0481 


1 


.52446- 


•005 


3 


.0532 


128 


2, 


.93540- 


■006 


3, 


.0267 


1, 


.8692c- 


-006 


3 


.0278 


256 


3, 


.6341e- 


•007 


3, 


.0139 


2, 


.31376- 


•007 


3 


.0141 


512 


4, 


.5205e- 


•008 


3, 


.0071 


2 


.87796- 


•008 


3 


.0071 


32 


1, 


.6061e- 


•005 






1, 


.01916- 


•005 






64 


9, 


.3375e- 


•007 


4, 


.1044 


5 


.94986- 


•007 


4 


.0983 


128 


5, 


.61640- 


■008 


4, 


.0553 


3 


.57506- 


-008 


4, 


.0568 


256 


3, 


.4266e- 


•009 


4, 


.0348 


2, 


.18156- 


•009 


4, 


.0346 


512 


1, 


.9952e- 


•010 


4, 


.1022 


1 


.27026- 


•010 


4, 


.1022 
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shows the space-time evolution of the exact quantity Qs and the temporal evolution of q 
totally averaged over Vis with different layering approaches; it can be seen that again 
the CVA approach seems to overestimate the advection speed, the same situation 
happens with the DGPO scheme, being the FEM source term slightly better than the 
PS on generating a scheme that tries to recover Qs] such performance is confirmed 
by Figures I4.3II4.41 where the ADER4-FEM scheme is the one exhibiting the lower 
residuals. 

4.2. 2D atmospheric model. We also study the performance of our approach 
in a linear, non-hydrostatic atmospheric model. If we start from conservation princi- 
ples, the model is usually expressed in terms of the velocity field in the x — z plane 
{u,w), the density p, the pressure p and the temperature T. A non-conservative, 
atmospheric-related formulation can be cast in terms of the potential temperature 9 
and the Exncr function tt which are defined as 

, n^(P-Y^^\ (4.6) 
\PoJ 

where Cp is the specific heat of dry air at constant pressure, R is the gas constant for 
dry air and po is atmospheric pressure at ground level. We can linearize this model 
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upon a decomposition into reference states and perturbations of the form ^{x, z, t) = 
^{z) + ^'{x,z,t), and assuming a reference isothermal state ( T constant), we can 
obtain a rescaled model for the perturbations 



dtq' + Ad^q' + B d,q' + Sq' ^0, 



where 



A 



u 




















' \ 





u 
















-M 


-L 








u 










M 








Cs 








u ) 




I 


c 





/ 



(4.7) 



(4.8) 



and 



B = 



/ ^ 

c, 



V c, y 



(4.9) 



We refer the reader to [12] for the complete derivation of this set of equations. The 
unknowns of the presented model are the velocity field perturbation in the x-z plane 
{u',w') , the potential temperature perturbation 9' and the Exncr function perturba- 
tion tt', all of them represented in the unknown vector q' — {u' , w' , 9', tt')*. Parameters 
in the isothermal atmosphere are usually initialized as follows. We prescribe a con- 
stant value for the reference temperature T and the Brunt- Vaisala frequency Af. The 
reference state for the potential temperature is then given by the relation 




9{0) = T. 



(4.10) 



The reference state for the Exner function is obtained via the hydrostatic balance 
relation 



-9, 



7f(0) = 1. 



(4.11) 
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Fig. 4.3. Linear advection and averaging: vertically averaged solutions (right) and its residuals 
j ^exact _ ^approx | ^^^^^^ ^ elements. Top: ADER4-FEM. Bottom: ADER4-PS 
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The sound speed Cg is constant, and can be computed from 



-Rife, 



(4.12) 



with Cy specific heat of dry air at constant volume, and the reference state for the 
density is obtained from the state equation of ideal gas 



R9 



(4.13) 



The remaining parameter U corresponds to the reference state for the horizontal 
velocity. For all the tests in this section, we set T ~ 300[K], N' = 0.01[s~^], g = 
9.8[ms-% Cy = 717[J kg~^ K-'^l Cp = 1006[J kg-^ R-^], R = 287 M[J kg~'^ R-^], 
Pq = 10^[Pa\ and U = 20[ms-^]. 
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Fig. 4.4. Linear advection and averaging: vertically averaged solutions (right) and its residuals 
^^exact _ ^approx ^ 100 x 20 elements. Top: ADER^-CVA. Bottom: ADER4-DGP0 
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So far we have presented a method for the resolution of scalar equations; with 
this example we extend it to constant-coefficient linear hyperbolic systems. As the 
dimensionally reduced model only treats it as an unidimensional system of balance 
laws in the x-direction, it is convenient at this point to transform the system into 
characteristic variables in that direction. Expressing the equations in terms of the 
characteristic vector q leads to 



(4.14) 



where 
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A 



( U + c, 



V 










\ 



U 

Q U J 



/ ic, \ 
ic, 


V c, y 

(4.15) 



and 



S =Q-^SQ = 



( 







\ -C -C -N J 



i(7r' + u') 



(4.16) 



For this example we consider three different vertical discretizations: piecewise linear 
continuous finite elements, and piecewise constant and hnear discontinuous Galerkin 
expansions. We close the system with periodic boundary conditions at a; = L^^ for 
every variable, and we set w' — at z — and z — Hz- 

In the case of the piecewise linear continuous Galerkin method, the extension to 
the system is straightforward; after uniformly partitioning the vertical interval into 
N elements, every variable is expanded as a linear combination of the same basis 
functions, and testing every equation separately leads to an assembled representation 
of the form: 









dtqh + A dx qii ^ 


where 








( 




\ 


A = 







{U-Cs)l 









Ul 




[ 





UI J 



, s = 



= 0, 






(4.17) 











CI \ 











CI 











2m 


2CI 


-2CI 


-2MI 


/ 








(4.18) 



and 



c = 

2 





































2M"iK 


2M-1K 









Qh = 



4 



1% 



(4.19) 



with I the identity matrix of size N x N, and M, K arc defined as in eqns. (|2.14p . 

For piecewise linear discontinuous Galerkin method (we omit the piecewise con- 
stant case as it is trivially contained) the only additional consideration is related with 
the definition of the upwind fiux for a linear system. For the i — th element the local 
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equations read: 



dt co,i + AdxCo,^ +S co.i + - {B q{z = Zi) - B q{z = Zi^i) } = 0, 



(4.20) 



dt ci.r + A 9, ci,, + 5 ci,, + -{B q{z = z,)+B q{z = - 2B cq^, } = 0. (4.21) 

In this case, being B a constant-coefficient matrix that admits a diagonal represen- 
tation Af, such that B = R^^AtR, we spHt A?, into A{|" and Aj^, containing only 

non-negative and non-positive eigenvalues respectively. We define B^ = R^^A^R 

and B = R^^Aj^R, and therefore the upwinding is defined by 



6g(z-z,_i) = 6 g,_i(C = l)+6 g.(e = -l), 



(4.22) 



6 9(z = z,_i) = S g,(e = l) + S g,+i(e = -l), 



(4.23) 



with 



B 



/-I 



4 



V 



-1 
-1 


2 



1 \ 

1 



-2 J 



B 



2i 

4 



/I 

1 


V 2 



1 
1 




1 





2/ 



(4.24) 



4.2.1. Convergence rates. Our first test is performed in order to validate the 
theoretical convergence rates of the ADERm scheme for hyperbolic balance laws in 
the presence of a stable source term. For this purpose we set an initial potential 
temperature perturbation 



9'{x,z,0) 



lOn 



X sm 



27r 



(4.25) 



and every other perturbation is initialized at rest; values of other test parameters are 
Lx = 300[fcm] and = 10[A:m]. We use a CFL number CFL — 0.5, and the time 
step is limited by the fastest wave of the model, i.e., 



At = CFL 



Aa 



U 



(4.26) 



being all the simulations stable in this case. Observation time is at t = 500 [s]. 
We consider a fixed number of 20 vertical elements. Results for relative errors and 
convergence rates in space at final time of simulation, with respect to a reference 
solution with 1024 horizontal elements, for FEM, DGPO and DGPl discretizations 
are shown in Tables l¥!^4. 41 We note that, excepting for the scheme with source term 
arising from a DGPO discretization, convergence rates in Li and Loo norm, for spatial 
order 2, 3 and 4 are in accordance to the expected values, i.e., the presence of the 
source terms arising from vertical discretizations is not a drawback to recover the 
same orders as in the purely advective problem. Henceforth further numerical studies 
are performed only with optimally convergent schemes, i.e., with FEM and DGPl 
source terms. 
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Table 4.3 

Convergence rates for the ADER scheme with a FEM source term for the atmospheric model. 



Method 


N 


Loo error 


order 


Li error 


Li 

i 


order 


ADER2-FEM 


32 


1.2854e-002 




1.3350e-002 








04 


z.oioue-uuz 




z.zUoze-uuz 








izo 


4.yuuze-uuo 




4.oUD4e-UUo 




z, 


.zoyo 




ZOD 


i.ooo ( e-UUo 


i.O ( (4 


i.irjDoe-uuo 


1, 


.yo ( z 






6. ( oyoe-uu4 


1 CI KA 


z.yUo4e-UU4 


1, 


.yo /y 


ADER3-FEM 


32 


3.0623e-002 




3.05486-002 








64 


3.5159e-003 


3.1227 


3.5144e-003 


3, 


.1197 




128 


3.5367e-004 


3.3134 


3.53636-004 


3 


.3130 




256 


3.3499e-005 


3.4002 


3.35026-005 


3 


.3999 




512 


3.4205e-006 


3.2919 


3.42046-006 


3, 


.2920 


ADER4-FEM 


32 


2.9463e-003 




2.93016-003 








64 


4.0365e-004 


2.8677 


4.04166-004 


2, 


.8580 




128 


4.3203e-005 


3.2239 


4.31966-005 


3, 


.2259 




256 


3.0715e-006 


3.8141 


3.07196-006 


3, 


.8137 




512 


1.8900e-007 


4.0225 


1.88996-007 


4, 


.0227 



4.2.2. Inertia-gravity waves in a periodic channel. Having verified the 
numerical accuracy of the scheme, we turn our attention to a test case which is 
meaningful from an atmospheric modelling point of view. We study the behavior of 
non- hydrostatic scale gravity waves (NHIGW). In the same framework of 19 and 114) . 
with the same parameters as in the previous section, we test with an initial potential 
temperature perturbation given by: 



'ix,z,0) = 0.01- 



1 



TTZ \ 



(4.27) 



Computational domain and other parameters are given in Table 14.51 This test 
is performed in order to verify proper speed propagation of the inertia-gravity waves 
and its symmetry with respect to the advected initial perturbation. The initial per- 
turbation is expected to be translated by the background horizontal velocity, while at 
the same time breaking and radiating symmetrically to the left and to the right; the 
axis of symmetry should recover the background translation speed. 
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Table 4.4 

Convergence rates for the ADER scheme with a DGPO and DGPl source terms for the atmo- 
eric model. 



Method N error L^o order Li error Li order 



ADER2-DGP0 
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.2264e- 


■004 


1, 


.5549 


2 
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-004 


1 


.5542 


32 


1, 


.3685e- 


■002 






6 


.6063e- 


-003 






64 


6, 


.5137e-003 


1, 


.0710 


1 


.9915e- 


■003 


1 


.7300 
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Table 4.5 

Computational domain and parameters for simulations of NHIG W test case 



Test a Tf CFL Aa; h At 

NHIGW 300[/cm] 10[fcm] 5[fcm] 3 x 10^[s] 0.5 l[km] 500[m] lM[s] 

As in the convergence rates validation, we first observe that for every method 
a reduction of the CFL number to .5 is needed in order to have stable simulations. 
We conjecture that this is due to the presence of a source term which includes an 
oscillatory behavior that is reflected in a set of purely imaginary eigenvalues on its 
spectrum, independently of the nature of the vertical discretization. Although nonlin- 
ear, the treatment that the ADER scheme gives to the source term can be interpreted 
as a single-step, explicit method and therefore a reduction of the time step is required 
in order to include a section the imaginary axis inside the stability region of the 
scheme. Table shows values for the potential temperature extrema that are in ac- 
cordance with the results obtained for the NHIGW test case in [14] . Simulations show 
that the phenomena is accurately described in physical terms as it preserves symme- 
try in the propagation of the initial perturbation and the correct advection speed for 
both FEM and DGPl source terms (see l4.6lC71 and l?!5)) . These second-order vertical 
schemes manage to conserve the quantity u'^ + w'^ as shown in figure 14.51 which is 
not only of physical meaning, but also illustrates the stability of the scheme; finally, 
the velocity field and its divergence (see 14.81 and 14. 9p are qualitatively correct in the 
sense that they reflect oscillations for the divergence and "fan-like" behavior for the 
velocity field. 

Table 4.6 

Extrema att = 3, 000[s] for the NHIGW test case with ADER4 scheme for horizontal resolution, 
FEM and DGPl schemes for vertical discretization 



Source 0' 0' ?/' 7/' m' ?// tt' tt' 

^max min ^max min ^raax ^min max min 

FEM 2.70e-3 -1.43e-3 2.66e-3 -2.66e-3 7.69e-4 -8.93e-4 5.28e-4 -5.28e-4 

DGPl 2.63e-3 -1.36e-3 2.59e-3 -2.58e-3 7.18e-4 -8.38e-4 4.54e-4 -4.33e-4 



5. Summary and outlook. We have presented an approach for the generation 
of layered, high horizontal resolution advective models, which has been tested in 
linear advection models. We have implemented different approaches for the layering 
procedure, while the resulting dimensionally reduced system is solved via the ADER 
finite volume scheme with prescribed order of accuracy. In the test cases, the method 
has shown a good performance in terms of stability, convergence, and recovery of 
the underlying physics of the equations. A feature of the method is efficiency; the 
C-K procedure is computed only once before initialization, and the all the structure 
of the code can be easily parallelized/ vectorized. Moreover, the layered structure of 
the model can be parallelized in the context of domain decomposition. Thus, the 
extension of the presented procedure to very higher-order approximations seems to 
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Fig. 4.5. 1. Initial potential temperature perturbation for the NHIGW test case (top) 2. Evo- 
lution of u'^ + ui'^ during the NHIGW test case (bottom). 




be feasible and meaningful. Although we have limited our analysis to linear model, 
the discontinuous Galerkin layering approach can be easily extended to non-linear 
models due to its local character. The class of models where the presented approach 
can be applied is quite wide, being the main requirement that the advective behavior 
of the model is governed by the horizontal direction, which is a plausible assumption 
in most atmospheric phenomena. 
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